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We use recent theoretical advances to develop a new functional form for interatomic forces in 
bulk silicon. The theoretical results underlying the model include a novel analysis of elastic prop- 
erties for the diamond and graphitic structures and inversions of ab initio cohesive energy curves. 
The interaction model includes two-body and three-body terms which depend on the local atomic 
environment through an effective coordination number. This formulation is able to capture success- 
fully: (i) the energetics and elastic properties of the ground state diamond lattice; (ii) the covalent 
re-hybridization of undercoordinated atoms; (iii) and a smooth transition to metallic bonding for 
overcoordinated atoms. Because the essential features of chemical bonding in the bulk are built 
into the functional form, this model promises to be useful for describing interatomic forces in silicon 
bulk phases and defects. Although this functional form is remarkably realistic by usual standards, 
it contains a small number of fitting parameters and requires computational effort comparable to 
the most efficient existing models. In a companion paper, a complete parameterization of the model 
is given, and excellent performance for condensed phases and bulk defects is demonstrated. 
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I. INTRODUCTION 

The study of materials properties is increasingly rely- 
ing on a microscopic description of the underlying atomic 
structure and dynamics. While many of the key features 
can be described by a small number of atoms that are ac- 
tively participating in a physical process, many problems 
of interest require of order lO'^-lO^ or even higher num- 
ber of atoms and time scales of 10-100 ps for a proper de- 
scriptian. Ah initio methods based on density functional 
theoryEl and the local density approximation (DFT/LDA) 
have been intensively and successfully used tOjjprovide 
a microscopic description of simple structuread. For 
more complex cases, including for instance disordered or 
stepped surfaces, dislocations, grain boundaries, crystal 
growth and the amorphous-to-crystal transition, a large 
number of atoms is required, making an ah initio de- 
scription untenable. A possible alternative for these cases 
might be empirical interatomic potentials which are com- 
putationally much less expensive. The difficulty in em- 
ploying empirical potentials is their unproven ability to 
capture the physics of structures far from the fitting data 
used to construct them. Developing reliable empirical po- 
tentials remains an issue of great interest and possibly of 
great rewards. 

Silicon is a test case for the development of empirical 
potentials for covalent materials. Its great technologi- 
cal importance, the vast amount of relevant experimental 
and theoretical studies available, and its intrinsic interest 
as the representative covalent material make it an ideal 
candidate for exploring to what extent the empirical po- 
tential approach can be exploited. In recent years, more 
than 30 empirical potentials for silicon have been devel- 
oped and applied to a number of different systems, and 



more recently compared to each othem'H. They differ in 
degree of sophistication, functional form, fitting strategy 
and range of interaction, and each can accurately model 
various special atomic configurations. Sijifaces and small 
clusters are the most difficult to handlcfl'El, but even bulk 
material (crystalline and amorphous phases, solid defects 
and the liquid phase) has resisted a transferable descrip- 
tion by a single potential. Realistic simulations of impor- 
tant bulk phenomena such as plastic deformation, diffu- 
sion and crystallization are still problematic. 

In this article, we derive a general model for the func- 
tional form of interatomic forces in bulk tetrahedral semi- 
conductors. This functional form is appliedrio the proto- 
typical case of silicon in a companion articlaS. The devel- 
opment of the model is organized as follows: In section 0, 
we briefly review existing potentials and approximations 
of quantum models for silicon and extract important con- 
clusions about the desirable features of a successful in- 
teratomic potential. Recent theoretical advances used in 
deriving our model from ah initio total energy data are 
outlined in section III. A functional form that incorpo- 



rates the theoretical results using a minimal number of 
fltting parameters is presented and discussed in section 
[V. Finally, section ^ contains some concluding remarks. 



II. 



REVIEW OF EMPIRICAL POTENTIALS AND 
APPROXIMATIONS 

A. Empirical Potentials 



The usual approach for deriving empirical potentials is 
to guess a functional form, motivated by physical intu- 
ition, and then to adjust parameters to fit ah initio total 
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energy data for various atomic structures. A covalent 
material presents a difficult challenge because complex 
quantum-mechanical effects such as chemical bond for- 
mation and rupture, hybridization, metalization, charge 
transfer and bond bending must be described by an ef- 
fective interaction between atoms in which the electronic 
degrees of freedom have somehow been "integrated out"Ll. 
In the case of Si, the abundance of potentials in the liter- 
ature illustrates the difficulty of the problem and lack 
of specific theoretical guidance. In spite of the wide 
range of functional forms and fitting strategies, all pro- 
posed models possess comparable (and insufficient) over- 
all accuracyB. It has proven almost impossible to at- 
tribute the successes or failures of a potential to specific 
features of its functional form. Nevertheless, much can 
be learned from past experience, and it is clear that a 
well-chosen functional form is more useful than elaborate 
fitting strategies. 

To appreciate this point we compare and contrast some 
representative potentials for silicon. The pioneering po- 
tential of Stillinger and Weber (SW) has only eight pa- 
rameters and was fitted to a few experimental properties 
of solid cubic diamond and liquid siliconB. The model 
takes the form of a third order cluster potentialQ in which 
the total energy of an atomic configuration {Rij} is ex- 
pressed as a linear combination of two- and three-body 
terms, 

E = J2V2iR^J)+Y.V3{R,J,R,k), (1) 

ij ijk 

where Rij = Rj — Ri, Rij — \Rij \ and we use the conven- 
tion that multiple summation is over all permutations of 
distinct indices. The range of the SW potential is just 
short of the second neighbor distance in the equilibrium 
DC lattice, so the pair interaction V2{r) has a deep well 
at the first neighbor distance to represent the restoring 
force against stretching sp^ hybrid covalent bonds. The 
three-body interaction is expressed as a separable prod- 
uct of radial functions g(r) and an angular function h{9) 

Vz{rur2)^9[ri)g{r2)h{W2). (2) 

where I12 — cos6'i2 = fi ■ r2l(r\r2)- The angular func- 
tion, h(V) = (I + 1/3)^ , has a minimum of zero at the 
tetrahedral angle to represent the angular preference of 
sp? bonds, and the radial function g{j) decreases with 
distance to reduce this effect when bonds are stretched. 
The SW three-body term captures the directed nature of 
covalent sp^ bonds in a simple way that selects the dia- 
mond lattice over close-packed structures. Although the 
various terms lose their physical significance for distor- 
tions of the diamond lattice large enough to destroy sp^ 
hybridization, the SW potential seems to give a reason- 
able description of many states experimentally relevant, 
such as point defects, certaia surface structures, and the 
liquid and amorphous statesn. The SW potential contin- 
ues to be a favorite choice in the literature, due in large 



part to its appealing simplicity and apparent physical 
content. 

Another popular and innovative empirical model is 
th&.Terspff potential- with three versions generally called 
Tli, T2E2I, and Tsllil. The original version Tl has only 
six adjustable parameters, fitted to a small database of 
bulk polytypes. Subsequent versions involve seven more 
parameters to improve elastic properties. The Tersoff 
functional form is fundamentally different from the SW 
form in that the strength of individual bonds is affected 
by the presence of surrounding atoms. Using Carls- 
son's terminology, the Tersoff potential is a third or- 
der cluster functionalQ with the cluster sums appearing 
in nonlinear, pambinations. As suggested by theoretical 
argumentaiSnij, the energy is the sum of a repulsive pair 
interaction (^B.{r) and an attractive interaction 
that depends on the local bonding environment, which is 
characterized by a scalar quantity C, 

^ = Y. ['^flX^) + V{C,^i)'^A{R^,)\ (3) 

=E^3(%,i?^fc), (4) 

where the function represents the Pauling bond or- 
der. The three-body interaction has the form of Eq. (|^) 
with the important difference that the angular function, 
although still positive, may not have a minimum at the 
tetrahedral angle. The Tl, T2 and T3 angular functions 
are qualitatively different, possessing minima at 180°, 
90° and 126.745°, respectively. The Tersoff format has 
greater theoretical justification away from the diamond 
lattice than SW, but the three versions do not outperform 
the SW potential, overall, perhaps due to their handhng 
of angular forcestl. Nevertheless, the Tersoff potential is 
another example of a successful potential for bulk prop- 
erties with a physically motivated functional form and 
simple fitting strategy. 

The majcpitw-pf empiricai^otentials fall into either the 
generic SWll3il3 or Tersoflll3l23 formats just described, 
but there are notable exceptions that provide further in- 
sight into successful approaches for designing potentials. 
First, a number of potentials possess functional forms 
that have either limited validity or no physical motivation 
at all, suggesting that fitting without theoretical guid- 
ance risrjiot the optimal approaGlw-, The Valence Force 
Fieldc^lHJ and related poteatialsEScS (of which there are 
over 40 in the literatureEs) involve scalar products of 
the vectors connecting atomic positions, an approxima- 
tion that is strictly valid only for small departures from 
equilibrium. Thus, extending these models to highly 
distorted bonding environments undermines theip-ithe- 
oretical basis. The potential of Pearson ei. aZO, as 
the authors emphasize, is not physically motivated, but 
rather results from an exercise in fitting. Their use of 
Lennard-Jones two-body terms and Axilrod- Teller three- 
body terms, characteristic of Van der Waals forces, has no 
justification for covalent materials. The potential of Mis- 
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triotis, Flytzanis and Farantos (MFF)@ is an interesting 
attempt to include four-body interactions. Although the 
importance of four-body terms is certainly worth explor- 
ing, the inclusion of a four-body term in a linear cluster 
expansion is not unique, and-theoretical analysis tends 
to favor nonlinear functionala3'E£lEJ. 

A natural strategy to improve on the SW and Ter- 
soff models is to replace simple functional forms with 
more flexible ones and complement them with more elab- 
orate fitting schemes. The Holding and Andersen (BA) 
potentials generalizes the TersofF format with over 30 
adjustable parameters fit to an unusually wide range of 
structures. Although it has not been thoroughly tested, 
the BA potential appears to describe simultaneously bulk 
phases, defects, surfaces and .small clusters, a claim that 
no other potential can makeEI. However, its complexity 
makes it difficult to interpret physically, and since a large 
fitting database was used, it is unclear whether the po- 
tential can reliably describe structures to which it was 
not explicitly fit. In this vein, tiie spline- fitted potentials 
of the Force Matching MethodEj represent the opposite 
extreme of the SW and Tersoff approaches: physical mo- 
tivation is bypassed in favor of elaborate fitting. These 
potentials involve complex combinations of cubic splines, 
which have effectively hundreds of adjustable parame- 
ters, and the strategy of matching forces on all atoms in 
various defect structures is the most elaborate attempted 
thus far. Although the method may be worth pursuing 
as an alternative, it has not yet produced competitive 
potentialstiJ. Moreover, even if a reliable potential could 
result from such fitting strategies, it would make it hard 
to interpret the results of atomistic simulations in terms 
of simple principles of chemical bonding. Such interpre- 
tation is essential, in our view, if physical insight is to be 
gained from computer simulations. 

In spite of relentless efforts, no potential has demon- 
strated a transferable description of silicon in all its 
formgj leading us to another important conclusion: it 
may be too ambitious to attempt a simultaneous fit of 
all of the important atomic structures (bulk crystalline, 
amorphous and liquid phases, surfaces, and clusters) 
since qualitatively different aspects of bonding are at 
work in different types of structures. Theory and gen- 
eral experience suggest that the main ingredient needed 
to differentiate between surface and bulk bonding pref- 
erences is a more sophisticated description of the local 
atomic environment. A notable example in this respect is 
the innovative Thermodynamic Interatomic Force Field 
(TIFF) potential of Chelikowsky et. alE3, which includes 
a quantity called the "dangling bond vector" that is a 
weighted average of the vectors pointing to the neighbors 
of an atom. For symmetric configurations characteris- 
tic of the ideal (or slightly distorted) bulk material, the 
dangling bond vector vanishes (or is exceedingly small). 
Conversely, a nonzero value of the dangling bond vector 
indicates an asymmetric distribution of neighbors. While 
the TIFF danghng bond vector description appears to be 
very useful for undercoordinated structures like surfaces 



and small clusters, in this work we restrict ourselves to 
bulk material and thus use a simpler, scalar environment 
description. Our goal is to obtain the best possible de- 
scription of condensed phases and defects with a simple, 
theoretically justified functional form. 



B. Approximation of Quantum Models 



An alternative to fitting guessed functional forms 
is to derive potentials by systematic approximation of 
quantum-mechanical models. So far, this approach has 
failed to produce superior potentials, but important con- 
nections between electronic structure and effective inter- 
atomic potentials have been revealed. Although attempts 
are being made to directly approximate Density Func- 
tional Theorytj, the most useful contributions involve ap- 
proximating various Tight Binding (TB) models, which 
can themselves-be derived as approximations of first prin- 
ciples theorieaSJ. These methods are based on low or- 
der moment approximations of the TB local density of 
states (LDOS), which is used to express the avemgaJaajod 
energy as the sum of occupied bonding stateaJ'OESLj. 
Pettifor has derived a many-body potential, similar in 
form to the Tecsoff potential, by approximation of the 
TB bond ordeiO. More recently, an angular dependence 
remarkably close to the T3 angular function has been de- 
rived for a bonding from the lowest opder two-site term 
in the Bond Order Potential expansiorCJ, but the analyt- 
ically derived function has a flat minimum around 130° 
and thus differs qualitatively with the Tl and T2 poten- 
tials. With hindsight, a simple physical principle explains 
these results: a a bond is most weakened (desaturated) 
by the presence of an another atom when the resulting 
angle is small {9 < 100°) because in such cases the atom 
lies near the bond axis, thus interfering with the a orbital 
where it is most concentrated. Working within the same 
framework of the TB LDOS, Carlsson and coworkers have 
derived potentials with the Generalized Embedded Atom 
Methodc2rE3. Harrison has arrived at a similar model by 
expanding the average band energy in the ratio of the 
width of the bonding band to the bond-antibond split- 
ting, the relevant small parameter in semiconductors^^. 
These potentials resemble the SW potential in its descrip- 
tion of angular forces with an additive three-body term, 
particularly for small distortions of the diamond lattice. 
The transition to metallic behavior in overcoordinated 
structures involves interbond interactions similar to the 
Tersoff and embedded atom potentials. 

Many-body potentials can be derived from quantum- 
mechanical models if we restrict our attention to impor- 
tant small sets of configurations. Using a basis-tpf sp^ 
hybrid orbitals in a TB model, Carlsson et. aZ.Ql2a have 
shown that a generalization of the SW format, in which 
Eq. (^ is replaced by a ifQrm similar to that used by 
Biswas and Hamann (BH)ta, 
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V3{ri,r2) = ^ 9,n{ri)g„Ar2) 



(5) 
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is valid in the vicinity of the equihbrium diamond lattice. 
In general, the fourth moment controls the essential band 
gap of a semiconductor, implying four-body interactions, 
but the separable, three-body SW/BH terms are a conse- 
quence of the open topology of the diamond lattice: the 
only four-atom hopping circuit between firgt neighbors is 
the self-retracing path i^j^i^k^iU. 

We can make analogous arguments for the graphitic 
lattice to draw conclusions about sp^ hybrid bonds. Ig- 
noring the weak, long-range interaction between hexag- 
onal planes, we can assume a TB basis of sp^ hybrid 
orbitals and follow Carlsson's derivation. Because the 
self-retracing path is also the only first neighbor hopping 
circuit in a graphitic plane, a cluster expansion with the 
generic BH three-body interaction is also valid for hexag- 
onal configurations, with the functions in Eqs. (]l|) and 
differing from their diamond sp^ counterparts, as de- 
scribed below. These calculations also suggest that a 
locally valid cluster expansion should acquire strong en- 
vironment dependence for large distortions from the ref- 
erence configurationQ. 

These studies provide theoretical evidence that the lin- 
ear three-body SW/BH format is appropriate near equi- 
librium structures, while the nonlinear many-body Ter- 
soff format describes general trends across different bulk 
structures. For the asymmetric configurations found in 
surfaces and small clusters, these theories also suggest 
that a more complicated environment dependence than 
Tersoff's is needetl-like the dangling bond vector of the 
TIFF potentialEJO. In conclusion, direct approxima- 
tion of quantum models can provide insight into the ori- 
gins of interatomic forces, but apparently cannot pro- 
duce improved potentials. The reason may be that the 
long chain of approximations connecting first principles 
and empirical theories is uncontrolled, in the sense that 
there is no small parameter which can provide an asymp- 
totic bound foe, the neglected terms for a wide range of 
configuration^. 



III. INVERSION OF AB INITIO ENERGY DATA 

There are very few hard facts concerning the nature 
of interatomic forces. Although there has been a pro- 
liferation of ab initio energy and force calculations for 
a wide range of atomic structures, it has proven diffi- 
cult to discover any concrete information regarding the 
functional form of interatomic potentials. With the ubiq- 
uitous fitting approach, it is never clear whether discrep- 
ancies with ab initio data result from ea incorrect func- 
tional form or simply suboptimal fitting^. Thus, in addi- 
tion to the practical problem of designing potentials, it 
is also difficult to build a simple conceptual framework 
within which to understand the complexities of chemi- 
cal bonding. In this section, we summarize our recent 



efforts to extract features of interatomic forces directly 
from ab initio total energy data. In order to investigate 
the global trends in bonding across bulk structures pre- 
dicted by quantum theories, we first per form inversions of 
ab initio cohesive energy curves in part III A , B y ana lyz- 
ing elastic properties of covalent solids in part IIIB , we 
then explore the cohesive forces in certain special (high 
symmetry) bonding states, which can be viewed as an 
inversion of ab initio energies restricted to selected im- 
portant configurations. 



A. Inversion of Cohesive Energy Curves 

We have recently shown that it is possible to derive ef- 
fective interatomic [Potentials for covalent solids directly 
from ab initio datacjifj. The inversion procedure gener- 
alizes the "ab mitio pair potential" of Carlsson, Gelatt 
and EhrenreichE^ to many-body interactions and|Jor ar- 
bitrary strains beyond uniform volume expansion^. For 
the case of silicon, this work provides first principles evi- 
dence in favor of the generic bond order form of the pair 
interaction. 



1/2 (r, Z)^Mr)+p{Z)Mr), 



(6) 



where 4>R{r) represents the short-range repulsion of 
atoms due to Pauli exclusion of their electrons, 
represents the attractive force of bond formation, and 
p{Z) is the bond order, which determines the strength 
of the attraction as a function of the atomic environ- 
ment, measured by the coordination Z. Tl, 
cal behavior of the bond order is as foUowsL 
The ideal coordination for Si is Zq = 4, due to its va- 
lence. As an atom becomes increasingly overcoordinated 
{Z > Zq), nearby bonds become more metallic, charac- 
terized by delocalized electrons. In terms of electronic 
structure, the LDOS for overcoordinated atoms can be 
reasonably well described by its scalar second moment. It 
is a well established result that the leading order behavior 
of the bond ordet.is-J5(Z) ^ Z^^^^ in the second moment 
approximationlll'tj'Ea. For Z < Zq on the other hand, a 
matrix second moment treatment predicts cuxpughly con- 
stant bond order (additive bond strengths )E3. For small 
coordinations higher moments are needed to incorporate 
important features of band shape characteristic of co- 
valpnt bnTiHincr primarily the formation of a gap in the 
LDOS|3oOl3. Thus, the bond order should depart from 
the divergent Z~^l'^ behavior at lower coordinations with 
a shoulder at the ideal coordination oi Z — Zq where the 
transition to metallic Z^^^^ dependence begins. 

Inversion of ab initio cohesive energy curves verifies 
that trends in chemical bonding across various bulk 
bonding arrangementgL_are indeed consistent with these 
theoretical predictionscil. Previously, the only evidence 
in support of the bond order formalism came from equi- 
librium bond lengths and energies for a small set of ideal 
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crystal structuresBil^'tSllSI. The inversion approach has 
revealed for the first time that the bond order decompo- 
sition expressed by Eq. (^) is actually valid for a wide 
range of volumes away from equilibrium and for a repre- 
sentative set of low energy crystal structures. In addition 
to selecting the generic form of the pair interaction, in- 
version provides a precise measure of the relative bond 
orders in various local atomic configurations. For ex- 
ample, the bond order of sp^ bonds involving three-fold 
coordinated atoms is about 5% greater than that of four- 
fold coordinated sp^ bonds in silicon. 

These results have immediate implications for em- 
pirical potentials. The main conclusion is that the 
generic Tersoff format is much more realistic than the 
SW format for highly distorted configurations. However, 
the inversion results also indicate that a coordination- 
dependent pair interaction can provide a fair description 
of high-symmetry crystal structures without requiring 
additional many-body interactions. In particular, angu- 
lar forces are only needed to stabilize these structures 
under symmetry- breaking distortions, primarily for small 
coordinations. In order to make a quantitative connec- 
tion between Tersoff's functional form and our inverted 
ah initio data, angular contributions to the bond order 
must somehow be suppressed for ideal crystal structures. 

The inversion procedure applied to explicit three-body 
interactions has also led to somCj-useful conclusions. Al- 
though it is not always the casec3, inverted three-body 
radial functions g{r) tend to be strictly decreasing func- 
tions (like SW), especiallv— jvhen an overdetermined set 
of input structures is usecO. Inverted angular functions 
h{l) also tend to penalize small angles {9 < tt/2) less 
than most existing models, in agraement with a compar- 
ative study of empirical potentialsO. We must emphasize, 
however, that the results of this section concern general 
trends in chemical bonding, and have little to offer in 
terms of the precise nature of interatomic forces in spe- 
cial atomic configurations, such as the low-energy states 
of hybrid covalent bonds. To understand better these 
critical cases, we employ a related inversion strategy. 



B. Analysis of Elastic Properties 

A useful theoretical approach to guide the develop- 
ment of pottentials, which has been pursued recently only 
by CowleycJ, is to predict elastic properties implied by 
generic functional forms and compare with experimen- 
tal or ab initio data. This tool for understanding inter- 
atomic forces dates back to the 19th century, when St. 
Venant showed that the assumption of central pairwise 
forces supported by Cauchy and Poisson implied a reduc- 
tion in tha number of independent elastic constants from 
21 to 15113. The corresponding six dependencies, given 
by the single equation C12 — C44 if atoms are at cen- 
ters of ciftbi|e-|Symmetry, are commonly called the Cauchy 
relational3'c3. They provide a simple test for selecting 



which materials can be described by a pair potentialSi. 
Once it was realized that the Cauchy relations are not 
satisfied by the experimental data for semicondjUCjtors, a 
number of authors in this century, led by BornEHlEil, de- 
rived generalized Cauchy. relations for nonccntral forces 
in the diamond structurd£3't3. Building upon this body 
of work, we have recently analyzed the elastic proper- 
ties of several general classes of many-body potentials 
in the diamond and graphitic crystal structures in order 
to gain insight into the mechanical behasaor of sp^ and 
sp^ hybrid covalent bonds, respectiveljcJ. These high 
symmetry atomic configurations must be accurately de- 
scribed by any realistic model of interatomic forces in a 
tetravalent solid. Here we will only outline results di- 
rectly related to the model presented in the next section. 

sp^ Hybrids: Consider one of the simplest many-body 
interaction models for a tetrahedral solid, that is the 
generic SW format defined in Eqs. (|^) and (||), with near- 
est neighbor interactions and an angular function having 
a minimum of zero at the tetrahedral angle {h — ^ 
0, h" > 0). In that case, first considered by HarrisoroEj, 
the functional form of the potential has only two degrees 
of freedom for elastic behavior, V2' and h", the curva- 
tures of the pair interactioH-.and of the angular function 
at their respective minimacj. Since cubic symmetry al- 
lows for three independent elastic moduli, there is an 
implied relation, due to Harrison, 

{7Cn + 2Ci2)C44 - 3(Cii + 2C,2){Cn - C12). (7) 

Using the experimental datall shown in Table I, the ratio 
of the two sides of the Harrison relation is 1.16, indicating 
a reasonable description by a simple SW model. In con- 
trast, the potenjtials with the Tersoff format, T2, T3 and 
Dodson (D0D)ll3, are far from satisfying this relation. 
This does not imply rejection of the Tersoff format, be- 
cause the functional form has more than enough degrees 
of freedom to exactly reproduce all the elastic constants. 
However, as such, the inability of Tersoff potentials to 
accurately describe elastic behavior when constrained to 
fit other important properties does suggest a potential 
shortcoming in the functional form. 

A more compelling reason to select the SW format over 
others in the literature comes from the unrelaxed shear 
modulus C44 which does not include relaxation of the 
internal degrees of freedom in the crystal unit cell. In 
the early literature on elastic forces, unrelaxed elastic 
moduli were ignored, because they are not experimen- 
tally accessible. With the advent of ab initio calculations 
that predict elastic constants within a few percent of ex- 
perimental values, we can now analyze unrelaxed elastic 
properties as well. Considering again the simple SW for- 
mat, with its two degrees of freedom, we have discovered 
another relation for the unrelaxed moduli. 



4Cii + 5Ci2 = 9C?4 



(8) 



As shown in Table I, the experimental and ab initio elas- 
tic moduli satisfy this relation within experimental and 
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computational error. On the other hand, more general 
cluster potentials and functionals, including the TersofF 
format, BH and PTHT, do not require this relation, and 
in fact cannot satisfy it under the usual circumstances. 
This is demonstrated in Table I and explains why it has 
proven difficult tp-pbtain good elastic properties with the 
Tersoff potentiaE3. These results unambiguously select 
the SW format with first neighbor interactions for de- 
scribing small homogeneous strains of the diamond lat- 
tice. Although imperfect, internal relaxation with the 
SW format is also much better than with other models. 
Combining Eqs. (|^) and (||), we arrive at a relation in- 
volving all four moduli, Cn, C12, C44 and C44, 



AA ^4 
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8Ci2)2 



9(7Cii + 2Ci2)' 



(9) 



that expresses the effect of internal relaxation. If the two 
degrees of freedom in the SW format are used to repro- 
duce the experimental values of C\\ and C12, and thus 
also C?4 by Eq. (^) , then the predicted value of C44 from 
Eq. d) is 0.71 Mbar, which is only 12% smaller than 
the experimental value of 0.81 Mbar. The elastic behav- 
ior of the SW format is quite remarkable considering it 
has only half of the necessary degrees of freedom, while 
most other models are overdetermined for elastic behav- 
ior. This explains the surprising facta that the SW poten- 
tial gives one of the best descriptions of elastic properties 
in spite of not having been fit to any elastic constants. 
We conclude that it is the superiority of the simple SW 
functional form that gives the desirable properties, not a 
complex fitting procedure. 

Using analytic expressions for the elastic constants it 
is possible to devise a simple prescription to achieve good 
elastic properties with the SW format. As a simple con- 
sequence of 1/3) = 0, the curvature of the pair po- 
tential is given by. 



^"(^'i) = TT(^ii+2Ci2). 



(10) 



The curvature of the angyilar function can be related to 
the second shear modulu£0. 



gM2;i"(-l/3) = ^(Cn-Ci2) 



(11) 



where r^, ad and Vrf — denote the equilibrium first 
neighbor distance, lattice constant and atomic volume. 
Using the ah initio data in Table I, the right hand sides 
of Eqs. ^ and (0) evaluate to 8.1 eV/l^ and 5.7 eV, 
respectively. This provides a simple two-step procedure 
to maintain good elastic behavior while fitting any model 
with the SW format near the diamond lattice: {%) scale 
the pair interaction V2 ir) to obtain the correct bulk mod- 
ulus K = (di -I- 2Ci2)/3, and (ii) scale the three-body 
energy to set the second shear modulus. As shown above, 
this will lead to perfect unrelaxed elastic constants and 
only a 12% error in C44. 



sp^ Hybrids: We have also obtained useful informa- 
tion about interatomic forces due to sp^ hybrid iionds 
from the elastic moduli of the graphitic structureEll. In 
this analysis we neglect interplanar interactions, which 
are insignificant compared to the covalent bonds within 
a single, hexagonal plane. Our goal is to understand the 
elastic properties of sp^ hybrids appearing around three- 
fold coordinated atoms in a bulk enviroament, such as 
a dislocation core or a grain boundarjiS. An isolated 
hexagonal plane embedded in three-dimensional space 
has three independent (unrelaxed) elastic constaats, Cn, 
C12, and C44 with units of energy per unit areaE^. It can 
be shown that C44 = for any three-body cluster poten- 
tial or functional, in perfect agreement with the vanishing 
ab initio valueEJ. There is no relation for the remaining 
constants, Cn and C12, implied by empirical models be- 
cause each functional form possesses at least two degrees 
of freedom. 

Drawing on the TB approximations described above, 
which correctly predict the general form of interactions 
mediated by sp^ hybrids, we proceed by assuming a sep- 
arate three-body cluster potential for sp^ hybrids given 
by Eqs. (||) and (||). By analogy with the sp'^ case, 
we further assume the simpler SW form of Eq. (||) for 
the three-body interaction, with the important difference 
that the angular function has a minimum of zero at the 
hexagonal angle of 27r/3 rather than at the tetrahedral 
angle. We again restrict the interaction range to nearest 
neighbors engaged in the covalent bonds that dominate 
cohesion. These are not the only possible choices, but 
we can evaluate their validity through analysis of elastic 
moduli. _ 

With such a functional foHaEil, which differs from all 
existing empirical potentialgl£3, stability considerations 
imply Cn > 3Ci2, which is indeed satisfied by the Mh 
initio values, Cn = 1-79 Mbar and C12 = 0.51 MbaiH. 
More importantly, we can relate the mechanical proper- 
ties of sp^ and sp^ hybrids. The relative radial stiffness 
is given by a simple ratio of elastic constants. 



9rl Vd{Cn + 2Cu) 



(12) 



where the subscript h refers to the equilibrium hexag- 
onal plane with area per atom = a^\/3/4, and d 
refers to the diamond lattice. Using the ah initio re- 
sult, rh = 2.23A, the prefactor in Eq. (|l|) is 0.99, so 
the elastic constant ratio in parentheses provides a direct 
comparison of sp^ and sp'^ radial forces. The ab initio 
value of that ratio is 1.4 ± 0.1, implying that sp^ bonds 
have 40% greater radial stiffness than sp'^ bonds. The 
same result also follows directly from inverted oair po- 
tentials for the graphitic and diamond structurestd. 

A similar elastic analysis yields an expression for the 
relative angular stiffness of sp^ and sp'^ hybrid bonds, 



fe;:(-l/2) ^ 256gd{rd)^ AhjCn - 3Ci2)h 
h'^{-l/3) 243gh{rh)^ Ki(Cn - Cis)^ ' 



(13) 
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Using the ab initio data, we have the general result, 
ghirhfK{~l/2)/gdirdYh'^{-l/3) = 0.46 ±0.15. As- 
suming gd{r) ~ ghi^) with each_£nnction decreasing in 
accordance with inversion resultalj, then the product of 
prefactors in Eq. (^) is nearly unity. In that case the ra- 
tio of elastic constants in parentheses allows us to quan- 
tify the relative bending strength of the hybrid bonds. 
The ab initio value for the ratio of 0.44 ± 0.15 indicates 
that the angular stiffness of sp^ bonds is smaller than 
that of sp^ bonds by about a factor of two, in spite of 
the greater radial stiffness of sp^ bonds. Our conclusion 
for the relative bending strength of sp^ and sp'^ hybrids 
would be reversed only \l gg(rg) were smaller than gdird) 
by at least a factor of two, which seems unlikely in light 
of the bond orders. 

Elastic constant analysis suggests that a hybrid co- 
valent bond is well represented by a separable, first- 
neighbor, three-body cluster potential whose angular 
function has a minimum of zero at the appropriate angle. 
This may seem to contradict the ample evidence we have 
cited in favor of the Tersoff format for large distortions of 
the diamond lattice, particularly those involving changes 
in coordination. These findings are consistent, however, 
in light of Carlsson's argument that cluster potentials like 
SW can accurately fit narrow ranges of configurations 
while cluster functionals like Tersoff's provide a less ac- 
curate but physicalhL. acceptable fit to a much broader 
set of configurationsEj. 

This body of results forms a reliable foundation upon 
which to build empirical potentials for bulk tetravalent 
solids. In general, we conclude that the functional form 
of atomic interactions should reduce exactly to appropri- 
ate cluster potentials in special bonding geometries, with 
environment dependence that interpolates smoothly be- 
tween these special cases and captures general trends. 
We shall refer to this theoretically motivated functional 
form as the Environment Dependent Interatomic Poten- 
tial (EDIP) for Bulk Si. 



IV. FUNCTIONAL FORM 

Although reasonable interaction potentials can be de- 
rived using the analytic methods of the previous section, 
such inversion schemes become most powerful when used 
as theoretical guidance for fitting. The reason is that 
inversion necessarily involves a restricted set of ab ini- 
tio data. Although the input data can be perfectly re- 
produced (unless it is overdetermined), it is desirable to 
allow an imperfect description of the inversion data in 
order to achieve a better overall fit of a wider ab initio 
database that includes low symmetry defect structures. 
Thus, our approach is to incorporate the theoretically 
derived features of the previous section directly into our 
functional form, and then to fit the potential to a care- 
fully chosen ab initio database with a minimal number 
of parameters. In this way, we can systematically derive 



a reliable potential for bulk properties while keeping the 
functional form simple enough to allow for efficient com- 
putation of forces as well as intuitive understanding of 
chemical bonding in covalcnt solids. 



A. Scalar Environment Description 

The simplest description of the local environment of 
an atom is the number of nearest neighbors, determined 
by an effective coordination number Zi for atom i, 

Z, - /(^™) (14) 

where f{Ri,n) is a cutoff function that measures the con- 
tribution of neighbor m to the coordination of i in terms 
of the bond length Rim ■ The special sp^ and sp^ bonding 
geometries can be uniquely specified by their coordina- 
tions due to their high symmetry. Since environment 
dependence is not needed in those cases, it is natural to 
take the coordination number to be a constant, except 
when large distortions from equilibrium occur. Moreover, 
covalent bonds tend to involve only first neighbors, as in- 
dicated by ab initio charge densityi-ealculations of open 
structures like the diamond latticeE3. Thus, we choose 
the neighbor function to be exactly unity for typical co- 
valent bond lengths, r < c, with a gentle drop to zero 
above a cutoff b that excludes second neighbors, 

{1 if r < c 

exp(^3^) if c<r<6 . (15) 
if r > 6 

where x — (r c)/(b ~ c). This particular choice of cut- 
off function is appealing because it has two continuous 
derivatives at the inner cutoff c, and is perfectly smooth 
at the outer cutoff b. The cutoffs b and c are restricted to 
lie between first and second neighbors of both the hexag- 
onal plane and diamond lattice in equilibrium, so that 
their coordinations are 3 and 4, respectively. 

Our scalar description of the atomic environment is 
similar to Tersoff's, but there are notable differences. 
First, the perspective is that of the atom rather than 
the bond: With our potential, the preferences for spe- 
cial bond angles, bond strengths and angular forces are 
the same for all bonds involving a p|ar| ti c ju| p; atom. This 
is in contrast to the Tersoff formatErllilll3'E3 in which a 
mixed bond-atom perspective is adopted: the contribu- 
tion of atom i to the strength of bond (ij) is affected by 
the "interference" of other bonds (ik) involving atom i. 
This model provides an intuitive explanation for trends 
in chemical reaction paths of moleculesEil and allows for 
both covalent and metallic bonds to be centered at the 
same atom, as observed, for example, in ab initio charge 
densities for the BCT5 latticeE3, which lies between the 
covalent diamond lattice and the metallic /3-tin lattice. 
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However, the analysis of elastic properties discussed ear- 
lier favors the present approach for environment depen- 
dence near the diamond lattice. Another important dif- 
ference between our model and Tersoff's is the separa- 
tion of angular dependence from the bond order. As 
we shall see, this allows us to control independently the 
preferences for bond strengths, bond angles, and angu- 
lar forces in a way that the Tersoff potential cannot. By 
keeping the bond order simple, we can also directly use 
the important theoretical results that motivated the Ter- 
soff potential in the first place. 



B. Coordination-Dependent Chemical Bonding 

Our potential consists of coordination-dependent two- 
and three-body interactions corresponding to the defin- 
ing features of covalent materials: pair bonding and an- 
gular forces. The energy of a configuration is a sum 
over single-atom energies, E = Ei, each expressed as 
a sum of pair and three-body interactions 



E^ =J2V2{R^J,Z^)+Y^ Vs {R^o , R^k , Z, 



(16) 



jk 



depending on the coordination Zi of the central 
atom. The pair functional V2{Rij, Zi) represents the 
strength of bond (ij), while the three-body functional 
V'i{Rij,Rik,Zi) represents preferences for special bond 
angles, due to hybridization, as well as the angular forces 
that resist bending away from those angles. From our 
atomic perspective, the pair interaction is broken into a 
sum of contributions from each atom, and similarly the 
three-body interaction is broken into a sum over the three 
angles in each triangle of atoms. Note that due to the 
environment dependence, the contributions to the bond 
strength from each pair of atoms are not symmetric in 
general, ^2(^,^1) 7^ V2{Rji,Zj). 

Pair Bonding: We adopt the well-established bond or- 
der format of Eq. (H) for the pair interaction. Drawing 
on the popularity of the SW potential, we use those func- 
tional forms for the attractive and repulsive interactions. 



V2{r,Z) = A 



cxp 



(17) 



which go to zero at the cutoff r = a with all deriva- 
tives continuous. This choice can rcwoduce the shapes 
of inverted pair potentials for siliconcil. Because we have 
constructed Z, and hence p{Z), to be constant near the 
diamond lattice, our pair interaction reduces exactly to 
the SW form for configurations near equilibrium, thus 
allowing us to obtain excellent elastic properties as ex- 
plained above. Making this choice of repulsive term with 
the parameters obtained by fitting to defect structuresEI, 
we can follow the procedure of Bazant and KaxirasEJ to 
extract the implied bond order p{Z) from ab initio co- 
hesive energy curves for the following crystal structures 



(with coordinations given in parentheses): graphitic (3), 
diamond (4), BC-8 (4), BCT-5 (5), /3-tin (6), SC (6) and 
BCC (8). These structures span the full range from three 
and four-fold coordinated covalent bonding in sp^ and 
sp'^ arrangements, to overcoordinated atoms in metallic 
phases. The inverted ab initio bond order versus coordi- 
nation is shown in Fig. 1, along with two additional data 
points. Since we have only first neighbor interactions in 
the diamond lattice, we can obtain another bond order 
for three-fold coordination from the ab initio formation 
energy (3.3 eV) for an unrelaxed vacancy. An additional 
data point for unit coordination comes from the experi- 
mental binding energy|-(3.24 eV) and bond length (2.246 
A) of the Si2 moleculeEa. 

The bond order data has a clear shoulder at Z = Zq = 
4 where the predicted transition from covalent to metallic 
bonding occurs. For overcoordinated atoms with Z > Zq, 
the bond order approaches its rough asymptotic behav- 
ior, p oc characteristic of metallic band struc- 
ture. For coordinations Z < Zq, the bond order departs 
from the divergence, due to the formation of a 
band gap in the LDOS associated with covalent bonds. 
A natural choice to capture this shape is a Gaussian, 
p{Z) = e~^^ . In Fig. 1, we see that the bond or- 
der function we obtain from fittingEl is fairly close to the 
inversion data. It is intentionally somewhat too large 
for coordinations 5-8 to compensate for the small, but 
nonvanishing many-body energy for those structures, as 
described below. The collapse of the attractive functions 
0^(r) = (V2(r, Z)-VAir))/p{Z) with this choice of bond 
order shown in Fig. 2 is reasonably good, thus justifying 
the bond order formalism across a wide range of volumes. 
Our potential is the first to have a bond order in such 
close agreement with theory, which is a direct result of 
our novel treatment of angular forces. 

Angular Terms: In a thorough comparative study of Si 
potentials, Balamane et. al. attribute the limitations of 
empirical models to the inadequate description of angular 
forcescl. Our potential contains a number of innovations 
in handling angular forces, leading to a significant im- 
provement over existing models in reproducing ab initio 
data. Analysis of elastic properties shows that, at least 
near equilibrium, the three-body functional should be ex- 
pressed as a single, separable product of a radial function 
g{r) for both bonds and an angular function h{0, Z), 



V^{Rij,Rik,Zi) — g{Rij)g{Rik)h{lijk, Zi 



(18) 



Although the radial functions could vary with coordina- 
tion, in the interest of simplicity we have focused on the 
angular function as the most important source of coordi- 
nation-dependence. Inversion of ab initio cohesive energy 
curvealJ suggests that a consistent choice for the radial 
functions is the monotonic SW form, 



exp 
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-b 



(19) 
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which also goes to zero smoothly at a cutoff distance b, 
a value different from the two body cutoff a. Having 
separate cutoffs for two and three-body interactions is 
reasonable because they describe fundamentally differ- 
ent features of bonding. Although the pair interaction 
might extend considerably beyond the equilibrium first 
neighbor distance, the angular forces should not be al- 
lowed to extend beyond first neighbors, if they are to be 
interpreted as representing the resistance to bending of 
covalent bonds. 

Much of the new physics contained in our potential 
comes from the angular function h{l,Z). Theoretical 
considerations lead us to postulate the following general 
form: 

where H{x), w{Z) and t{Z) are generic functions whose 
essential properties we now describe. The overall shapc-of 
the angular function is given by H{x), a nonnegativeES'Ll 
function with a quadratic minimum of zero at the origin, 
H{0) = H'{0) = and H"{0) > 0. The function H{x) 
should also become flat away from the minimum well at 
the origin, resulting in zero angular force for large distor- 
tions away from equilibrium. This feature, which is ab- 
sent in most potentials including SW, is essential for the 
angular term to have physical meaning far from equilib- 
rium. When covalent bonds are strongly bent from their 
equilibrium angle they are weakened and replaced by new 
electronic states. Thus, for large angular distortions it is 
not possible to define a restoring force that drives atoms 
back towards an equilibrium bond angle. These proper- 
ties of H(x) can be satisfied by the following choice, 

i/(.x)=A(l-c--'), (21) 

which is similar in shape to the MFF angular functionS. 
However, our angular dependence is considerably more 
sophisticated than MFF due to its environment depen- 
dence. 

Motivated by theory, we choose the function t{Z) to 
control the coordination-dependent minimum of the an- 
gular function^JftlZ) = cos{6o{Z)) = —t{Z), with the 
following fornj£^'E3, 

t{Z) ^ui+ U2{u^e-^^^ - e-""^^). (22) 

The parameters, ui = —0.165799 , U2 — 32.557, uz — 
0.286198, and U4 — 0.66, were chosen to make the pre- 
ferred angle 9o{Z) — cos~^[— r(Z)] interpolate smoothly 
between several theoretically motivated values, as shown 
in Fig. 3: We have already argued that r(4) — 1/3 
and r(3) = 1/2 (so that sp^ and sp^ bonding corre- 
spond to the diamond and graphitic structures respec- 
tively), which determines two of the four parameters in 
r(z). The remaining two parameters were selected so 
that t(2) = t(6) = or 6lo(2) = 6io(6) = 7r/2. For 



two-fold coordination, this choice reproduces the prefer- 
ence for bonding along two orthogonal p-states with the 
low energy, nonbonding s state fully occupied. For six- 
fold coordination, the choice 6*0(6) — 7t/2 also reflects 
the p character of the bonds. However, structures with 
Z — 6 like SC and /3-tin are metallic, with delocalized 
electrons that tend to invalidate the concept of bond- 
bending underlying the angular function, a crucial point 
we shall address shortly. The vanishing many-body en- 
ergies for the graphitic plane and diamond structures al- 
low fitting of the pair interactions V2(r, 3) and V2(r, 4) 
to be guided by Eq. ([lO|), which determines V2"(rd,4) 
from the bulk modulus, and Eq. (p2|), which requires 
^2"(^/i>3)/F2"(rrf,4) w 1.4. Moreover, the shifting of the 
minimum of the angular function in our model incorpo- 
rates coordination-dependent hybridization in a way that 
other potentials cannot. 

Through the function w{Z), our angular function has 
another novel coordination dependence to represent the 
covalent to metallic transition. The width of the min- 
imum w(Z) is broadened with increasing coordination, 
thus reducing the angular stiffness of the bonds as they 
become more metallic. Similarly, as coordination is de- 
creased from 4 to 3, the width of the minimum is in- 
creased to reproduce the smaller angular stiffness of sp^ 
bonds compared to that of sp^ bonds. Thus, the function 
'w{Z) should have a minimum at Zq = 4 and diverge with 
increasing Z . Fitting of the model can be guided by Eq. 
dill), which determines w{A) from the second shear mod- 
ulus, and by Eq.(|l3|), which requires w{3)/w{A) « \/2. 
The softening of the angular function is important be- 
cause it allows the decrease in cohesive energy per atom 
concomitant with overcoordination to be modeled by a 
weakening of pair interactions. In contrast, cluster po- 
tentials like SW penalize overcoordinated structures with 
increased three-body energy that overcomes the decrease 
in pair bonding energy. This is an unphysical feature, 
since overcoordinated structures do not even have cova- 
lent bonds, and the many-body energy cannot be viewed 
as a consequence of stretching sp^ bonds far from the 
tetrahedral geometry. In this sense, the reasonably good 
description of liquid Si (a metal with about 6 neighbors 
per atom) with the SW potential appears to be fortu- 
itous. 

The coordination dependence of our angular function 
makes it possible for the first time to reproduce the well- 
known behavior of the bond order. The reason is that the 
contribution of the three-body functional to the total en- 
ergy is suppressed for ideal crystals and overcoordinated 
structures. The shifting of the minimum makes the three- 
body energy vanish identically for sp^ and sp'^ hybrids, 
and the variable width greatly reduces the three-body 
energy in metallic structures. With the three-body en- 
ergy suppressed, we can use our knowledge of the bond 
order for the graphitic, diamond, /3-tin and other lattices 
from inversion of cohesive energy curves to capture the 
energetics of these structures in the pair interaction, as 
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described above. Several other potentials have tried to 
incorporate the bond order predicted from theory, but 
the uncontrolled many-body energy makes it impossible 
to connect directly with theory. Our treatment of an- 
gular forces is intuitively appealing because the forces 
primarily model the bending of covalcnt bonds, with the 
control of global energetics left to the pair interactions. 

Although our model contains a complicated environ- 
ment dependence, forces can still be evaluated with com- 
putational speed comparable to much simpler existing 
potentials. The coordination dependence introduces an 
extra loop into each force calculation. For the three- 
body functional, this introduces a fourth nested loop over 
atoms m outside each triplet {ijk) that contribute to co- 
ordination of atoms i , which would make force evaluation 
much slower than the typical three-body cluster expan- 
sions used in most other potentials. However, our choice 
of /(r) greatly reduces the frequency of four body com- 
putations because nonzero forces result if the fourth atom 
lies in the range of being a partial neighbor, c < rim < 
which happens only for a small number of neighbors in 
most cases. If coordinations stay relatively constant dur- 
ing a simulation, as in a low temperature solid, the four- 
body force computation is insignificant. Indeed, we have 
found that force evaluation with jaur model can be almost 
as fast as with the SW potentialQ, which is an advantage 
of our model over others of comparable sophistication. 



sp^ and sp^ hybridization and some metallic states, but 
might also include more general situations in which atoms 
are more or less symmetrically distributed, like the liq- 
uid and amorphous phases and reconstructed dislocation 
cores and grain boundaries. The theory behind the model 
begins to break down for noninteger coordinations, since 
our effective coordination number is a way of smoothly 
interpolating between well-understood local structures. 
More seriously, no attempt is made to handle asymmet- 
ric distributions of neighbors, which are abundant in sur- 
faces and small clusters. Theory suggests that our model 
may be fitted to provide a good description of condensed 
phases and defects in bulk tetrahedral semiconductors, 
such as Si, Ge and with minor extensions perhaps alloys 
such as SiGe, that can be understood in terms of simple 
principles of covalent bonding. 
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V. CONCLUSION 



In summary, we have used recent theoretical innova- 
tions to arrive at a functional form that describes the de- 
pendence of chemical bonding on the local coordination 
number. Bond order, hybridization, metalization and an- 
gular stiffness are all described in qualitative agreement 
with theory. Consistent with our motivation, we have 
kept the form as simple as possible, reproducing the es- 
sential physics with little more complexity than exist- 
ing potentials. The fitted implenientation of the model 
described in the companion papeiO involves only 13 ad- 
justable parameters. Using the results of the present arti- 
cle, we provide theoretical estimates of almost half of the 
parameters, thus greatly narrowing the region of param- 
eter space to be explored during fitting. The remaining 
parameters are chosen to fit important bulk defect struc- 
tures. 

Considering the theory behind our model, we can an- 
ticipate its range of applicability. We have shown that 
the structure and energetics of the diamond lattice can be 
almost perfectly reproduced. Because small distortions 
of sp^ hybrids are accurately modeled, we would also 
expect a good description of the amorphous phase. De- 
fect structures involving sp^ hybridization should also be 
well described. In general, the model should perform best 
whenever the coordination number can adequately spec- 
ify the local atomic environment. This certainly includes 
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TABLE I. Comparison of elastic constants (in units of Mbar) for diamond cubic silicon computed from eipeirical models with 
experimental or ab initio (LDA) values. The values for experiment (EXPT) are from Simmons and WangLJ, for tight-binding 
(TB) from Bernstein and KaxirasEJ and for the empirical potentials Biswas-Haman (BH), Tersoff (T2, T3), Dodson (DOD) 
and Pearson- Takai-Haliciogl|tt-|Tiller (PTHT) from Balamanaj. The Stillinger-Weber (SW) values were calculated with the 
analytic formulae of Cowleyllj and scaled to set the binding energy to 4.63 eVQ. In the lower half of the table, we test the 
elastic constant relations discussed in the text by calculating the ratios an = (7Cii + 2Ci2)C44/3(Cii + 2Ci2)(Cii — C12) and 





EXPT 


LDA 


SW 


BH 


T2 


T3 


DOD 


PTHT 


TB 


Cii 
C12 
C44 
O44 


1.67 
0.65 
0.81 


1.11 


1.617 
0.816 
0.603 
1.172 


2.042 
1.517 
0.451 
1.049 


1.217 
0.858 
0.103 
0.923 


1.425 
0.754 
0.690 
1.188 


1.206 
0.722 
0.659 
3.475 


2.969 
2.697 
0.446 
2.190 


1.45 
0.845 
0.534 
1.35 


Oh 

Cub 


1.16 

0.99 




1.00 
1.00 


0.98 
1.67 


2.99 
1.10 


2.31 
0.89 


1.69 
0.27 


1.71 
1.29 


2.80 
0.82 



13 



1.2 



1.1 



0.9 



0.8 



0.7 



0.6 



1 


■, 


1 1 




1 






DIA 








VAC 


BC8 \ 










'o. 

BCT5 


pTIN ■.. 








1 1 




BCC ^ 

1 



1 



FIG. 1. Ab initio values for the bond order as a function of coordination, obtained from the inversion of cohesive energy 
curves for the graphitic (GRA), cubic diamond (DIA), BC8, BCT5, SG, /3-tin and EGG bulk structures and with additional 
points for the unrelaxed vacancy (VAC) and the dimer molecule (Si2), as explained in the text. For comparison the solid 
line shows the Gaussian p{Z) obtained from fitting to defect structures. The dotted line shows the dependence, the 

theoretically predicted approximate behavior for large coordinations. 
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FIG. 2. Attractive pair interactions from inversion of ab initio cohesive energy curves for the structures in Fig. 1 using 
the bond order and repulsive pair potential of our model. The solid lines are for the covalent structures with coordinations 
3 and 4, while the dotted lines axe for the overcoordinated metallic structures. The reasonable collapse of the attractive pair 
potentials indicates the validity of the bond order functional form of the pair interaction across a wide range of volumes and 
crystal structures. 



14 




FIG. 3. The coordination dependence of the preferred bond angle 9o{Z) (in degrees), which interpolates the theoretically 
motivated points for Z = 2,3, 4, 6, indicated by diamonds. 
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